Computing phase diagrams for a quasicrystal-forming patchy-particle system 
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We introduce an approach to computing the free energy of quasicrystals, which we use to calculate phase dia- 
grams for systems of two-dimensional patchy particles with five regularly arranged patches that have previously 
been shown to form dodecagonal quasicrystals. We find that the quasicrystal is a thermodynamic ally stable phase 
for a wide range of conditions and remains a robust feature of the system as the potential's parameters are varied. 
We also demonstrate that the quasicrystal is entropically stabilised over its crystalline approximants. 



Quasicrystals are a type of aperiodic crystal structure: they 
have well-ordered structures, but cannot be described by a pe- 
riodic lattice, not even one with incommensurate lattice param- 
eters [1]. They were first reported by Shechtman et al. [2], in 
whose work a rapidly cooled Al-Mn alloy was found to result in 
a metastable quasicrystal. Most quasicrystals discovered so far 
are metallic alloys [1]; however, recently, an increasing number 
of examples have been reported in the field of soft condensed 
matter [3]. Quasicrystals have also been observed to form in 
computer simulations, in which not only binary, but also one- 
component quasicrystals have been seen [4-8]. While certain 
quasicrystals have been shown to be stable in experiment [9], 
in simulations, the relatively short timescales accessible mean 
that it is not necessarily clear whether quasicrystals are truly 
the stable phase, or rather a metastable phase that is more ki- 
netically accessible. In this work, we address this important 
question of assessing the stability of quasicrystals and apply 
it to a soft matter system that has the potential to be realised 
experimentally. 

To prove that quasicrystals are thermodynamically stable, 
their free energy must be computed. However, such a calcula- 
tion is not straightforward: the principal problem is that there 
is no obvious reference state whose free energy is known and 
from which thermodynamic integration [10, 11] could be used 
to calculate the free energy of the quasicrystalline phase. A 
phase transition intervenes when integrating from an ideal gas, 
which is normally used as a reference state for fluid phases, 
and an Einstein crystal, used as a reference state for crystalline 
phases, would also not be suitable because it would fail to 
capture the configurational entropy associated with the qua- 
sicrystal's many possible structures. 

One approach to assessing quasicrystal stability is to com- 
pute the free energy of an approximant crystalline phase. How- 
ever, it is not immediately obvious whether quasicrystals are 
more or less stable than their approximants [5, 6]; whilst the 
enthalpy and vibrational properties of the quasicrystal and its 
approximant phases are likely to be similar, the free energy of 
a quasicrystal has a contribution from its configurational en- 
tropy that is not present for the approximant [7]. Nevertheless, 
such calculations have provided some insight into the phase 
behaviour of systems such as hard tetrahedra and bipyramids 
[6] and spherical micelles [7]. By contrast, Kiselev et al. re- 
cently estimated the free energy of the quasicrystal itself by 
combining a phonon contribution for a particular quasicrystal 
configuration from thermodynamic integration and a configura- 



tional contribution based on an approximation of uncorrelated 
phason flips [12]. They confirmed that, within their approxi- 
mation, the quasicrystal's free energy is in certain conditions 
lower than the approximant's. 

Our solution to this conundrum is to note that thermody- 
namic integration is not the only way the melting point of a 
solid can be computed. Another method is to simulate the 
direct coexistence of two phases with an interface [10]. By 
performing such simulations at a range of temperatures, we 
can bracket the regions where the quasicrystal melts (T > 7f us ) 
and grows (T < 7f us ). At T = 7f us , the chemical potential of 
the quasicrystal is equal to that of the fluid, which we can cal- 
culate using thermodynamic integration. Such an equilibrium 
approach by construction accounts for the quasicrystal's config- 
urational entropy. Once the quasicrystal's free energy is known 
at one point, we can use thermodynamic integration to reach 
other points of interest on the phase diagram. 

Here, we apply this approach to a two-dimensional system 
of patchy particles [13] with five regularly arranged attractive 
patches, which we model using a simple potential that has pre- 
viously been used in simulations of self-assembly and crystalli- 
sation [8, 14-16]. The potential is based on the Lennard-Jones 
(LJ) form, V L3 ( rij )=4e[(a L1 /nj) 2n - (a u /r y ) n ], where n = 6 
for the standard LJ potential, where for rn > Gu, this potential 
is multiplied by an angular modulation factor 
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where Op W is the 'patch width' (in radians) and Q^j is the angle 
between the patch vector of patch k (chosen to minimise this 
angle) on particle i and the interparticle vector 

This interaction potential provides a simple model for the 
patchy colloidal particles that many experimental groups have, 
with progressively increasing success, been seeking to develop 
[13, 17, 18]. The enhanced range of structural behaviour that 
such patchy interactions could facilitate has been extensively 
studied in computer simulation [8, 14-16, 19]. 

Dodecagonal quasicrystals were previously found to form 
on cooling systems of these five-patch particles in certain con- 
ditions, and structures with the lowest enthalpy were also iden- 
tified [8]. Some of the relevant phases studied are shown in 
Fig. l(a-d); particles are classified based on a common neigh- 
bour analysis [8] into a, H and Z environments (shown in 
Fig. 1(e)), by analogy to the Frank-Kasper phases [20]. The 
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FIG. 1. Examples of the main phases studied, (a) a phase, (b) Ap- 
proximate crystal, (c) Plastic hexagonal (Z) phase, (d) Quasicrystal, 
N = 2500. (e) Neighbour classification of a, H and Z environments. 
The numbers of common neighbours are given for each of the central 
particle's neighbours. These local environments are used to colour 
the particles in (a)-(d), where in addition, particles that do not adopt 
one of these environments are coloured green, (e) Diffraction pattern 
corresponding to the configuration shown in (a), (f) Dodecagonal 
motif characteristic of the quasicrystal. 

hexagonal (Z) phase (Fig. 1 (c)) was found to form at high pres- 
sures (as it is densest) and for wide patches (where the potential 
is closer to being isotropic). At low pressures and at reasonably 
narrow patch widths, the structure best satisfying the attractive 
patches is the enthalpically favoured phase. A crystal cannot 
exist with five-fold symmetry, and so no crystal with five per- 
fectly aligned patches exists. As a compromise, each particle 
has five neighbours, but the angles between the neighbours do 
not perfectly match the five-fold symmetry of the particle itself. 
Local environments satisfying this requirement are the o and H 
environments shown in Fig. 1(e); the lowest enthalpy structure 
at low pressure and reasonably narrow patch widths is the a 
crystal shown in Fig. 1(a) [8]. 

Quasicrystals were observed in some cooling runs in the re- 
gion where the a phase is the lowest in enthalpy, but near to the 
hexagonal 'boundary' . The structure of a typical quasicrystal 
is shown in Fig. 1(d); its diffraction pattern (Fig. 1(f)) exhibits 
a characteristic twelve-fold symmetry. The local environments 
in the quasicrystal are predominantly of the a type, but a signifi- 
cant fraction of Z environments can be seen: these are typically 



located at the centre of a dodecagonal motif (Fig. 1(g)). Much 
of the structure can be analysed in terms of packing of such do- 
decagons into triangular, square and rectangular arrangements 
[8], which can readily be seen in Fig. 1(d); the various possi- 
ble ways of arranging these dodecagons are likely to lead to 
substantial entropy. There is no translational order, but bonds 
can be orientated in any one of twelve directions, resulting in 
long-range orientational order of dodecagonal character. 

We perform Monte Carlo simulations in the NpT ensemble 
using a rectangular box with periodic boundaries [21]. As a 
starting point in the determination of phase diagrams for this 
system, we chose a patch width, temperature and pressure 
at which the quasicrystal is known to form on cooling, and 
thus locate the fluid-quasicrystal (F-QC) equilibrium transition. 
This transition is mostly rapid and facile, with essentially no 
hysteresis, which allows us to determine the coexistence points 
directly by performing simulations in relatively large boxes 
starting from both phases. Whilst we initially performed direct 
coexistence simulations to determine the melting point of the 
quasicrystal, these simulations mostly did not prove to be any 
more efficient than direct brute-force simulations, in which no 
initial interface was introduced, as the growth of the phases 
was not restricted to the initial interface. This lack of hysteresis 
is most likely indicative of a particularly low F-QC interfacial 
free energy. Phase transitions can be observed from a kink 
in the potential energy or density plotted against the variable 
we are changing (such as the temperature or the pressure). 
However, we can further test for the nature of the phase by 
calculating appropriate diffraction patterns: the quasicrystal 
exhibits a distinctive 12-fold diffraction pattern (Fig. 1(f)), the 
Z phase exhibits a six-fold pattern, and the fluid phase does not 
exhibit a well-defined pattern. 

Once a F-QC coexistence point is known, we can equate the 
chemical potential of the fluid (as calculated by thermodynamic 
integration) with that of the quasicrystal. From here, we can 
use thermodynamic integration to calculate the chemical poten- 
tial of the quasicrystal at other temperatures and pressures. By 
also performing thermodynamic integration to find the chemi- 
cal potential of the a phase, we can locate the putative a-QC 
transition. We plot in Fig. 2 the free energies of the o, QC 
and F phases as a function of T at a constant p/k^T . From 
this figure, we can identify two clear transitions (o-QC and 
QC-F), and a temperature window in which the quasicrystal is 
thermodynamically stable. In the same figure, we also show 
the free energy of an approximant crystal phase that is based on 
the most common packing of dodecagons in the quasicrystal; 
it is never thermodynamically stable. To rationalise why the 
quasicrystal is stable, we also calculate the enthalpy (obtained 
directly from simulations) and hence obtain the entropy as a 
function of the temperature. These functions are also depicted 
in Fig. 2, and it is clear that in the region in which the qua- 
sicrystal is stable, it does indeed have a higher enthalpy than 
the a phase, but its significantly higher entropy nevertheless 
allows it to become more stable. The origin of this entropy is 
a combination of the differing vibrational properties and the 
configurational disorder present in the quasicrystalline phase. 
If we assume the vibrational properties of the quasicrystal and 
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FIG. 2. Per-particle Gibbs energies, enthalpies and entropies, relative 
to the o" phase, as a function of T at o^p/k^T = 1.5, O pw = 0.49 
and n = 6. Lines end at the limit of (meta)stability of the phases. 
Dashed lines denote coexistence points, fcB7o«Qc/ e = 0.173 and 
^b7qc«f/ £ = 0.1995. The stable phase in each region is explicitly 
marked. The error in Ag/k^T is < 0.012^^ [21]. 

its approximant are identical, the configurational component 
of this entropy can be estimated from the entropy difference 
between the quasicrystal and its approximant, As/k^ « 0.113 
at r ooQC [22]. 

From Fig. 2, we can determine two coexistence points on 
the phase diagram. We construct the remainder of the phase 
diagram partly through the same procedure as above at, say, 
other pressures, but also through the use of Gibbs-Duhem inte- 
gration [10, 24]. This method allows us to integrate Clapeyron 
equation analogues to obtain new coexistence points from al- 
ready known coexistence points, but requires metastability and 
hysteresis in the transition so that the simulated phases retain 
their identity when simulated at and around the coexistence 
point. As a result, the method cannot be applied to the QC-F 
and QC-Z transitions for much of the phase diagram due to 
their relative reversibility. 

The resulting phase diagram is shown in Fig. 3(a). We note 
that the quasicrystal is stable for a wide range of pressures. 
At very high pressures, the Z phase becomes stable due to 
its greater density. The Z phase is a plastic crystal, meaning 
that particles are able to rotate relatively freely: this is under- 
standable, as there is no obvious preferred way to orientate a 
five-fold particle in a six-fold environment [25]. No orienta- 
tional ordering ensues in the Z phase even at low temperatures. 
Whether it would form an orientationally-ordered crystal or 
an orientational glass at sufficiently low T is not clear. It is 
noteworthy that for the QC-Z transition, dp/dT < 0; since 
Pz > Pqc, the Clapeyron equation implies that Sz > Sqc'- this 
is likely to stem from the orientational entropy of the plastic 
crystal. Interestingly, in addition to the two triple points (o- 
QC-Z and QC-Z-F), there is a range of pressures where the 
stable thermodynamic phase changes from F to Z to QC to a 
as the system is cooled. 
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FIG. 3. Phase diagrams. Markers show the technique used to obtain 
each point: Frenkel-Ladd and thermodynamic integration (filled cir- 
cles), Gibbs-Duhem integration (crosses), hamiltonian Gibbs-Duhem 
integration (triangles) and direct simulation (open circles). Points 
obtained from previous phase diagrams are shown as squares, (a) 
p/T-T phase diagram; C7 pw = 0.49, n = 6. (b) T-a pvj phase diagram; 
O^jp/ksT = 0.5, n = 6. (c) T-n phase diagram; o^p/k^T = 0.5, 
a DW = 0.49. 



We also look at the effect of modifying potential param- 
eters on the phase diagram, thus allowing us to assess the 
extent of quasicrystal stability and how finely-tuned the parti- 
cle properties would need to be to observe a quasicrystalline 
phase in experiment. To determine the phase diagram as the 
potential parameters change, we additionally use hamiltonian 
Gibbs-Duhem integration, whereby we numerically integrate 
a Clapeyron equation that has been generalised to allow for 
changes to the potential itself [10]. 

We first investigate the behaviour of the system as the patch 
width is varied. From the phase diagram in Fig. 3(b), we see 
that the quasicrystal is only stable for a limited range of a pw , 
although the range may perhaps be wider at other pressures. As 
the patch width narrows, the quasicrystal becomes increasingly 
enthalpically destabilised with respect to the a phase, as the 
many six-fold environments cannot satisfactorily fulfil patch- 
patch interactions. By contrast, as the patch width increases, 
the quasicrystal becomes more stable with respect to the o 
phase, but the hexagonal phase is stabilised even more: as 
the patches are wider and closer to the isotropic case, six-fold 
environments are enthalpically preferred. 

Another noteworthy feature of this phase diagram is that 
below Opw PS 0.35, r GO p begins to decrease: as the patches 
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become narrower, there is an increasing energetic penalty for 
patches that do not point directly at each other. For very nar- 
row patches, the a phase transforms into a distorted a phase 
(Da) [16], in which three of the five patches point directly at 
the patches of neighbouring particles, thus breaking the square 
symmetry of the lattice. The reason for this behaviour is that 
at narrow patch widths, only 'perfect' connections count, as a 
slight misalignment of the patches rapidly reduces the interpar- 
ticle attraction, and three perfect interactions become favoured 
over five imperfect ones. 

Because colloidal interactions are often quite short-ranged, 
we also wish to investigate what happens to the system when 
we change the LJ exponent («), which changes the range of 
the potential and the width of the potential well. We show the 
phase diagram as a function of n in Fig. 3(c). The effects of 
increasing n on isotropic potentials are well understood [26]: 
the liquid phase is energetically destabilised because the intrin- 
sic disorder in interparticle neighbour distances is penalised by 
the narrower potential [27]. In Fig. 3(c), the range of stability 
of the quasicrystal initially increases. The enthalpies of the QC 
and a phases increase only slightly as n increases, because their 
relative order means that most interparticle neighbour distances 
can be close to optimal, but do so more for the quasicrystal, as 
some disorder is typically present in the quasicrystal configu- 
rations (for example, particles with unclassified environments 
in Fig. 1(d)); the o-QC coexistence temperature does therefore 
increase slightly. By contrast, the fluid is initially much more 
enthalpically penalised: the fluid enthalpy increases rapidly, 
and the fluid density decreases as n increases. At the maximum 
in Tqcof (n ~ 14), the size of the temperature window of qua- 
sicrystal stability is three times that of the n = 6 LJ potential. 
However, a narrower potential well associated with increasing 
n also means that there is less vibrational entropy associated 
with the ordered structures, and beyond n 14, this effect be- 
comes dominant. The QC-F coexistence temperature decreases 
until the quasicrystal loses stability at n ~ 53. 

The potential at the higher values of n shown in Fig. 3(c) 
is very short-ranged indeed. The fact that the quasicrystal 
maintains and even increases its window of stability for ranges 
of attraction typical of colloidal particles suggests that this is 
a robust, general result and that stable quasicrystals can be 
expected to form in experimental realisations of this system. 

There are several potential approaches to such experimental 
realisations. If patchy colloidal particles similar to the ones 
studied here could be created, it would be reasonable to ex- 
pect that a quasicrystal may form when they are confined in 
two dimensions. For example, Chen et al. observed the for- 
mation of a 2D kagome lattice from tri-block Janus particles 
by introducing a density mismatch with the solvent to confine 
their colloidal system into two dimensions [18]. A possible 
alternative to using colloidal patchy particles might involve the 
use of DNA multi-arm motifs [28], for which a variety of 2D 
crystalline arrays have been observed for motifs with different 
numbers of arms, including a a phase for five-arm motifs that 
is analogous to what we see in the current model when patches 
are narrow. However, such DNA motifs have a well-defined 
'valence', and there is no equivalent of the patch width that 



could be tuned to make two co-ordination numbers compete. 
Therefore, a possible approach to forming a DNA quasicrys- 
tal might be to use a two-component mixture of five-arm and 
six-arm motifs of the appropriate composition. 

In summary, we have shown how we can compute if and 
where a quasicrystalline phase is thermodynamically stable. 
To the best of our knowledge, this is the first such calculation 
of the chemical potential of a quasicrystal and the associated 
phase diagrams that has been obtained directly from simula- 
tions with no approximations. For our patchy particle system, 
we found that the quasicrystalline phase is stable over a signifi- 
cant portion of the phase diagram and is stabilised primarily by 
its configurational entropy. It is robust to parameter changes 
in the model, which inspires confidence that such a thermody- 
namically stable quasicrystal might be experimentally realised. 

We wish to thank the Engineering and Physical Sciences 
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Supplementary material 



Simulation details 

The simulations described in our work used the Metropolis 
Monte Carlo approach [29] in the isobaric-isothermal ensemble 
[30, 31] with periodic boundary conditions. 

The potential as discussed in the manuscript is given more 
formally by 



V(r ih «Pi,<P;) = j ll^'L . C-) 



V L1 (r ij )V A (f lj ,( Pi ,( P j) OLj<r l7 , 



where r,j is the interparticle vector connecting the centres of 
the two particles i and j, rtj is its magnitude, <p, and (pj are the 
orientations of the particles i and j, the generalised Lennard- 
Jones potential is 

V L1 (r u ) = 4e [(a LJ /r, 7 ) 2 " - (ou/nj)"] , (3) 

where n = 6 for the standard Lennard- Jones potential, and 



V K {r ijl f h <P;)=exp 



-e 



2cJ p 2 w 



exp 



2 a 2 



(4) 



where Op W is the patch width, Q^j is the angle between the 
patch vector of patch k on particle i and the interparticle vector 
rn, and k m j n is the patch that minimises the magnitude of this 
angle. Two particles therefore interact only through a single 
pair of patches. In the simulations reported here, we use a 
potential cutoff of r cut = 3Clj, and the potential in Eq. (3) is 
shifted so that it equals zero at ru — r cut . The crossover to 
including angular modulation in the potential (Eq. (2)) is like- 
wise adjusted so that it still occurs when the potential energy 
is identically zero. 

Most simulations reported used relatively small system sizes 
(of the order of 500 particles), but the phase diagram in Fig. 3(a) 
and the free energies in Fig. 2 were also calculated using a 
system of the order of 2500 particles, and certain transitions 
were calculated with systems of 5000 and 10000 particles, 
with no detectable differences in the behaviour of the systems, 
confirming that reasonably small simulations were sufficient to 
calculate phase diagrams. 
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Summary of methods 



This can be integrated to give 



Here, we summarise the methods we used in the manuscript 
in somewhat more detail. The approaches described here are 
all standard and are provided here as a quick summary for the 
benefit of the reader only. 



Direct coexistence vs brute force simulation 

Direct coexistence simulations [10, 32] are simulations in 
which an explicit interface is set up between the two phases in 
question, and the thermodynamically stable phase is expected 
to grow into the thermodynamically unstable phase such that 
the interface between the two phases moves. In brute force sim- 
ulations, no interface is set up; instead, a single phase is simu- 
lated and it spontaneously converts into the thermodynamically 
stable phase. Typically, this brute force simulation process is 
prone to hysteresis effects, where for example the quasicrystal 
will form at some temperature T\ , but when the quasicrystal 
is heated in a reverse brute force simulation, it will melt at 
temperature T2, where T2 > T\, because nucleation involves a 
free energy barrier that must be overcome. In direct coexis- 
tence simulations, this hysteresis does not occur, because both 
phases are present from the start, and the interface between 
them is pre-formed. 

In our simulations, we established that there is essentially no 
hysteresis between the quasicrystal and the fluid under the con- 
ditions described in this paragraph, and so direct coexistence 
simulations offer no advantage over brute force simulations. 
Indeed, because nucleation is so facile, numerous phase inter- 
faces are spontaneously formed in a direct coexistence simula- 
tion, and so direct coexistence simulations were not useful in 
phase diagram determination at this stage. 

However, we have used direct coexistence simulations in 
certain cases as a check in confirming that the Frenkel-Ladd 
simulations yielded the correct coexistence points. 



r 2 



U 



G2 Gj 
Nk B T 2 ~ Nk B Ti J Tl Nk B T 2 



rdT, 



(7) 



and this is the result we require for thermodynamic integration. 
Provided we know the Gibbs energy at some T\, then we can 
find the Gibbs energy at T2 along the iso-(j3 p) curve by running 
several NpT simulations for temperatures between T\ and T2 
(at fixed /3 p), and finding a polynomial fit to the integrand, 
which can then by integrated analytically. Analogous results 
arise for integration along isobars [10], 



G 2 
Nk B T 2 

and isotherms [10], 

A 2 



H 



dT, 



Nk B T Nk B T 



Pl 



Pl 



k B Tp 



dp. 



(8) 



(9) 



It is also possible to extend the scheme by coupling to a 
different hamiltonian, such that the potential energy is given 
by U — U Ka i + XU ex tva, and then integrating along A; we find 
that [10] 



A (A = 1) =A(A=0)- 



(U e: 



dX. 



(10) 



This is known as hamiltonian thermodynamic integration [10] 
and is one component of the Einstein crystal Frenkel-Ladd 
approach [11]. 

The free energy of a fluid phase evaluated from the ideal gas 
by decreasing the pressure at constant temperature is given by 
[10] 



— -=ln(pA i )-l 
Nk B l Jo 



k B Tp z 



dp, 



(ID 



where A is the de Broglie thermal wavelength. It is usually 
beneficial to fit the entire integrand to a single polynomial 
to avoid diverging terms. As the density tends to zero, the 
integrand tends to the second virial coefficient [10]; it is helpful 
to evaluate this explicitly to benchmark the simulation results. 



Thermodynamic integration 



Thermodynamic integration [10, 11] refers to a series of inte- 
gration schemes of thermodynamic potentials. For integration 
along iso-(j3p) curves, where /3 = l/k B T, we start from the 
product rule 



(P» 



T \dT J 



(Pp) 



j2 ' 



(5) 



from where an application of simple thermodynamics allows 
us to write 



djG/T) 
dT 



(Pp) 



U 

f2- 



(6) 



Gibbs-Duhem integration 

In Gibbs-Duhem integration [10, 24], we numerically inte- 
grate the Clapeyron equation, 



dp_ 
dT 



AH 
TAV ' 



(12) 



or an analogous expression, in order to obtain new coexistence 
points from already known coexistence points. In our simu- 
lations, we used the fourth-order Runge-Kutta algorithm to 
perform this integration. 
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Frenkel-Ladd Einstein crystal approach 

The Einstein crystal approach [10, 11] is a fairly involved 
way of calculating the free energy of a crystal from the Einstein 
crystal, for which the free energy is known. We direct the 
interested reader to Frenkel and Smit's book (Ref. [30]) and 
the review paper by Vega and co-workers (Ref. [10]) for an in- 
depth discussion of the method. In the Frenkel-Ladd scheme 
as implemented for our model, we take the rotational energy to 
be 



= A m , sin 



{<Pi-<Pi,ang) 



(13) 



where A rot is a constant that we vary and p = 5 is the number 
of patches. We thus ensure that the particle's symmetry is 
accounted for: rotations into degenerate positions give the 
same Frenkel-Ladd energy. 



Reliability of free energy calculations 

Computing absolute free energies from simulations requires 
a combination of numerical techniques, each of which can 
introduce errors, both systematic and random. Since it is dif- 
ficult to keep track of how the errors propagate, the approach 
which is usually undertaken is to perform a posteriori checks 
on the whole process. As discussed by Vega and co-workers 
[10], 'consistency checks' can be applied to simulations to 
ensure that the data obtained in free energy calculations are 



self-consistent and correspond to reality. There are several 
such checks that can be performed. For example, we can cal- 
culate the free energy at a particular point using several routes 
(e.g. by integrating along different isotherms, isobars or iso- 
(/3 p) curves to the same point). Furthermore, direct coexistence 
simulations can be used to ensure that the transition tempera- 
ture does in fact occur where free energy curves intersect. Fi- 
nally, in simulations involving Gibbs-Duhem integration, we 
can integrate first in the 'forward' and then in the 'reverse' di- 
rection, ensuring that the method reproduces the starting point. 
We have performed all these checks to ensure that the results 
we have obtained are robust. 

Gao and co-workers compared the results of Debye crystal 
reference state calculations with different force constants to 
obtain an estimate of the error in the free energy of hexag- 
onal ice [33]. We do something similar to estimate the er- 
ror in free energy in Fig. 2 in the main paper: we calculate 
the Frenkel-Ladd free energy at several points (k B T /e = 0.1, 
k B T/e = 0.15, k B T/e = 0.17 and k B T/e = 0.19) along the 
iso-(j3/?) curve, and then use thermodynamic integration to 
obtain free energies from each of the starting Frenkel-Ladd 
free energies. The largest difference in free energy between 
any pair of calculations along the entire curve can then serve 
as an estimate of the error in free energy in our simulations. 
The largest such deviation was A(Ag)/k B T = 0.012, which is 
sufficiently small not to affect any conclusion, and the error in 
practice is likely to be significantly smaller, as we were able 
to confirm phase transition temperatures in independent direct 
coexistence simulations. 



